Skip to content

Conversation

bojle
Copy link
Contributor

@bojle bojle commented Aug 22, 2025

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.

@lntue @PiJoules

Copy link

Thank you for submitting a Pull Request (PR) to the LLVM Project!

This PR will be automatically labeled and the relevant teams will be notified.

If you wish to, you can add reviewers by using the "Reviewers" section on this page.

If this is not working for you, it is probably because you do not have write permissions for the repository. In which case you can instead tag reviewers by name in a comment by using @ followed by their GitHub username.

If you have received no comments on your PR for a week, you can request a review by "ping"ing the PR by adding a comment “Ping”. The common courtesy "ping" rate is once a week. Please remember that you are asking for valuable time from other developers.

If you have further questions, they may be answered by the LLVM GitHub User Guide.

You can also ask questions in a comment on this PR, on the LLVM Discord or on the forums.

@llvmbot
Copy link
Member

llvmbot commented Aug 22, 2025

@llvm/pr-subscribers-libc

Author: Shreeyash Pandey (bojle)

Changes

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.


Full diff: https://github.com/llvm/llvm-project/pull/154914.diff

6 Files Affected:

  • (modified) libc/config/linux/riscv/entrypoints.txt (+1)
  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/include/stdfix.yaml (+8)
  • (modified) libc/src/__support/fixed_point/fx_bits.h (+41)
  • (modified) libc/src/stdfix/CMakeLists.txt (+14)
  • (modified) libc/test/src/stdfix/CMakeLists.txt (+16)
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 1a03683d72e61..b783dd4b04c74 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -986,6 +986,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
     libc.src.stdfix.idivulr
     libc.src.stdfix.idivuk
     libc.src.stdfix.idivulk
+    libc.src.stdfix.rdivi
   )
 endif()
 
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5590f1a15ac57..f6beccc7229dd 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -1019,6 +1019,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
     libc.src.stdfix.idivulr
     libc.src.stdfix.idivuk
     libc.src.stdfix.idivulk
+    libc.src.stdfix.rdivi
   )
 endif()
 
diff --git a/libc/include/stdfix.yaml b/libc/include/stdfix.yaml
index 5b385124eb63d..451330c3478d2 100644
--- a/libc/include/stdfix.yaml
+++ b/libc/include/stdfix.yaml
@@ -544,3 +544,11 @@ functions:
     arguments:
       - type: unsigned long accum
     guard: LIBC_COMPILER_HAS_FIXED_POINT
+  - name: rdivi
+    standards:
+      - stdc_ext
+    return_type: fract
+    arguments:
+      - type: int
+      - type: int
+    guard: LIBC_COMPILER_HAS_FIXED_POINT
diff --git a/libc/src/__support/fixed_point/fx_bits.h b/libc/src/__support/fixed_point/fx_bits.h
index 00c6119b4f353..13eecba981664 100644
--- a/libc/src/__support/fixed_point/fx_bits.h
+++ b/libc/src/__support/fixed_point/fx_bits.h
@@ -224,6 +224,47 @@ idiv(T x, T y) {
   return static_cast<XType>(result);
 }
 
+inline long accum nrstep(long accum d, long accum x0) {
+  auto v = x0 * (2 - (d * x0));
+  return v;
+}
+
+/* Divide the two integers and return a fixed_point value
+ *
+ * For reference, see:
+ * https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
+ * https://stackoverflow.com/a/9231996
+ */
+template <typename XType> LIBC_INLINE constexpr XType divi(int n, int d) {
+  // If the value of the second operand of the / operator is zero, the
+  // behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16
+  LIBC_CRASH_ON_VALUE(d, 0);
+
+  unsigned int nv = static_cast<unsigned int>(n);
+  unsigned int dv = static_cast<unsigned int>(d);
+  unsigned int clz = cpp::countl_zero<unsigned int>(dv) - 1;
+  unsigned long int scaled_val = dv << clz;
+  /* Scale denominator to be in the range of [0.5,1] */
+  FXBits<long accum> d_scaled{scaled_val};
+  unsigned long int scaled_val_n = nv << clz;
+  /* Scale the numerator as much as the denominator to maintain correctness of
+   * the original equation
+   */
+  FXBits<long accum> n_scaled{scaled_val_n};
+  long accum n_scaled_val = n_scaled.get_val();
+  long accum d_scaled_val = d_scaled.get_val();
+  /* x0 = (48/17) - (32/17) * d_n */
+  long accum a = 2.8235lk; /* 48/17 */
+  long accum b = 1.8823lk; /* 32/17 */
+  long accum initial_approx = a - (b * d_scaled_val);
+  long accum val = nrstep(d_scaled_val, initial_approx);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  long accum res = n_scaled_val * val;
+  return static_cast<XType>(res);
+}
+
 } // namespace fixed_point
 } // namespace LIBC_NAMESPACE_DECL
 
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 843111e3f80b1..3cbabd17ad34c 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -89,6 +89,20 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
   )
 endforeach()
 
+foreach(suffix IN ITEMS r)
+  add_entrypoint_object(
+    ${suffix}divi
+    HDRS
+      ${suffix}divi.h
+    SRCS
+      ${suffix}divi.cpp
+    COMPILE_OPTIONS
+      ${libc_opt_high_flag}
+    DEPENDS
+      libc.src.__support.fixed_point.fx_bits
+  )
+endforeach()
+
 add_entrypoint_object(
   uhksqrtus
   HDRS
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index e2b4bc1805f7c..741522276feaa 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -120,6 +120,22 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
   )
 endforeach()
 
+foreach(suffix IN ITEMS r)
+  add_libc_test(
+    ${suffix}divi_test
+    SUITE
+      libc-stdfix-tests
+    HDRS
+      DivITest.h
+    SRCS
+      ${suffix}divi_test.cpp
+    DEPENDS
+      libc.src.stdfix.${suffix}divi
+      libc.src.__support.fixed_point.fx_bits
+      libc.hdr.signal_macros
+  )
+endforeach()
+
 add_libc_test(
   uhksqrtus_test
   SUITE

@llvmbot
Copy link
Member

llvmbot commented Aug 22, 2025

@llvm/pr-subscribers-backend-risc-v

Author: Shreeyash Pandey (bojle)

Changes

This PR includes only one of the fxdivi functions (rdivi). It uses a polynomial function for initial approximation followed by 4 newton-raphson iterations to calculate the reciprocal and finally multiplies the numerator with it to get the result.


Full diff: https://github.com/llvm/llvm-project/pull/154914.diff

6 Files Affected:

  • (modified) libc/config/linux/riscv/entrypoints.txt (+1)
  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/include/stdfix.yaml (+8)
  • (modified) libc/src/__support/fixed_point/fx_bits.h (+41)
  • (modified) libc/src/stdfix/CMakeLists.txt (+14)
  • (modified) libc/test/src/stdfix/CMakeLists.txt (+16)
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 1a03683d72e61..b783dd4b04c74 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -986,6 +986,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
     libc.src.stdfix.idivulr
     libc.src.stdfix.idivuk
     libc.src.stdfix.idivulk
+    libc.src.stdfix.rdivi
   )
 endif()
 
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5590f1a15ac57..f6beccc7229dd 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -1019,6 +1019,7 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
     libc.src.stdfix.idivulr
     libc.src.stdfix.idivuk
     libc.src.stdfix.idivulk
+    libc.src.stdfix.rdivi
   )
 endif()
 
diff --git a/libc/include/stdfix.yaml b/libc/include/stdfix.yaml
index 5b385124eb63d..451330c3478d2 100644
--- a/libc/include/stdfix.yaml
+++ b/libc/include/stdfix.yaml
@@ -544,3 +544,11 @@ functions:
     arguments:
       - type: unsigned long accum
     guard: LIBC_COMPILER_HAS_FIXED_POINT
+  - name: rdivi
+    standards:
+      - stdc_ext
+    return_type: fract
+    arguments:
+      - type: int
+      - type: int
+    guard: LIBC_COMPILER_HAS_FIXED_POINT
diff --git a/libc/src/__support/fixed_point/fx_bits.h b/libc/src/__support/fixed_point/fx_bits.h
index 00c6119b4f353..13eecba981664 100644
--- a/libc/src/__support/fixed_point/fx_bits.h
+++ b/libc/src/__support/fixed_point/fx_bits.h
@@ -224,6 +224,47 @@ idiv(T x, T y) {
   return static_cast<XType>(result);
 }
 
+inline long accum nrstep(long accum d, long accum x0) {
+  auto v = x0 * (2 - (d * x0));
+  return v;
+}
+
+/* Divide the two integers and return a fixed_point value
+ *
+ * For reference, see:
+ * https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
+ * https://stackoverflow.com/a/9231996
+ */
+template <typename XType> LIBC_INLINE constexpr XType divi(int n, int d) {
+  // If the value of the second operand of the / operator is zero, the
+  // behavior is undefined. Ref: ISO/IEC TR 18037:2008(E) p.g. 16
+  LIBC_CRASH_ON_VALUE(d, 0);
+
+  unsigned int nv = static_cast<unsigned int>(n);
+  unsigned int dv = static_cast<unsigned int>(d);
+  unsigned int clz = cpp::countl_zero<unsigned int>(dv) - 1;
+  unsigned long int scaled_val = dv << clz;
+  /* Scale denominator to be in the range of [0.5,1] */
+  FXBits<long accum> d_scaled{scaled_val};
+  unsigned long int scaled_val_n = nv << clz;
+  /* Scale the numerator as much as the denominator to maintain correctness of
+   * the original equation
+   */
+  FXBits<long accum> n_scaled{scaled_val_n};
+  long accum n_scaled_val = n_scaled.get_val();
+  long accum d_scaled_val = d_scaled.get_val();
+  /* x0 = (48/17) - (32/17) * d_n */
+  long accum a = 2.8235lk; /* 48/17 */
+  long accum b = 1.8823lk; /* 32/17 */
+  long accum initial_approx = a - (b * d_scaled_val);
+  long accum val = nrstep(d_scaled_val, initial_approx);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  val = nrstep(d_scaled_val, val);
+  long accum res = n_scaled_val * val;
+  return static_cast<XType>(res);
+}
+
 } // namespace fixed_point
 } // namespace LIBC_NAMESPACE_DECL
 
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 843111e3f80b1..3cbabd17ad34c 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -89,6 +89,20 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
   )
 endforeach()
 
+foreach(suffix IN ITEMS r)
+  add_entrypoint_object(
+    ${suffix}divi
+    HDRS
+      ${suffix}divi.h
+    SRCS
+      ${suffix}divi.cpp
+    COMPILE_OPTIONS
+      ${libc_opt_high_flag}
+    DEPENDS
+      libc.src.__support.fixed_point.fx_bits
+  )
+endforeach()
+
 add_entrypoint_object(
   uhksqrtus
   HDRS
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index e2b4bc1805f7c..741522276feaa 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -120,6 +120,22 @@ foreach(suffix IN ITEMS r lr k lk ur ulr uk ulk)
   )
 endforeach()
 
+foreach(suffix IN ITEMS r)
+  add_libc_test(
+    ${suffix}divi_test
+    SUITE
+      libc-stdfix-tests
+    HDRS
+      DivITest.h
+    SRCS
+      ${suffix}divi_test.cpp
+    DEPENDS
+      libc.src.stdfix.${suffix}divi
+      libc.src.__support.fixed_point.fx_bits
+      libc.hdr.signal_macros
+  )
+endforeach()
+
 add_libc_test(
   uhksqrtus_test
   SUITE

@PiJoules
Copy link
Contributor

Thanks. Could you add some tests for this?

@bojle
Copy link
Contributor Author

bojle commented Aug 22, 2025

Oops, forgot to add them. I've pushed the tests too. Please check now. @PiJoules

Copy link

github-actions bot commented Aug 26, 2025

⚠️ C/C++ code formatter, clang-format found issues in your code. ⚠️

You can test this locally with the following command:
git-clang-format --diff origin/main HEAD --extensions h,cpp -- libc/src/stdfix/rdivi.cpp libc/src/stdfix/rdivi.h libc/test/src/stdfix/DivITest.h libc/test/src/stdfix/rdivi_test.cpp libc/src/__support/fixed_point/fx_bits.h

⚠️
The reproduction instructions above might return results for more than one PR
in a stack if you are using a stacked PR workflow. You can limit the results by
changing origin/main to the base branch/commit you want to compare against.
⚠️

View the diff from clang-format here.
diff --git a/libc/src/stdfix/rdivi.cpp b/libc/src/stdfix/rdivi.cpp
index 227b41287..7021a8b1c 100644
--- a/libc/src/stdfix/rdivi.cpp
+++ b/libc/src/stdfix/rdivi.cpp
@@ -15,7 +15,7 @@
 namespace LIBC_NAMESPACE_DECL {
 
 LLVM_LIBC_FUNCTION(fract, rdivi, (int a, int b)) {
-  return fixed_point::divi<fract>(a,b);
+  return fixed_point::divi<fract>(a, b);
 }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/stdfix/DivITest.h b/libc/test/src/stdfix/DivITest.h
index f3ac145ad..8960c2cb5 100644
--- a/libc/test/src/stdfix/DivITest.h
+++ b/libc/test/src/stdfix/DivITest.h
@@ -44,14 +44,14 @@ public:
   }
 
   void testSpecial(DivIFunc func) {
-    EXPECT_EQ(func(0,10), 0.r);
-    EXPECT_EQ(func(0,-10), 0.r);
-    EXPECT_EQ(func(-32768,32768), FRACT_MIN);
-    EXPECT_EQ(func(32767,32768), FRACT_MAX);  
-    EXPECT_EQ(func(INT_MAX,INT_MAX), 1.0r);
-    EXPECT_EQ(func(INT_MAX-1,INT_MAX), 0.99999999r);
-    EXPECT_EQ(func(INT_MIN,INT_MAX), FRACT_MIN);
-    /* Expecting 0 here as fract is not precise enough to 
+    EXPECT_EQ(func(0, 10), 0.r);
+    EXPECT_EQ(func(0, -10), 0.r);
+    EXPECT_EQ(func(-32768, 32768), FRACT_MIN);
+    EXPECT_EQ(func(32767, 32768), FRACT_MAX);
+    EXPECT_EQ(func(INT_MAX, INT_MAX), 1.0r);
+    EXPECT_EQ(func(INT_MAX - 1, INT_MAX), 0.99999999r);
+    EXPECT_EQ(func(INT_MIN, INT_MAX), FRACT_MIN);
+    /* Expecting 0 here as fract is not precise enough to
      * handle 1/INT_MAX
      */
     EXPECT_EQ(func(1, INT_MAX), 0.r);
@@ -61,5 +61,5 @@ public:
 #define LIST_DIVI_TESTS(Name, XType, func)                                     \
   using LlvmLibc##Name##diviTest = DivITest<XType>;                            \
   TEST_F(LlvmLibc##Name##diviTest, Basic) { testBasic(&func); }                \
-  TEST_F(LlvmLibc##Name##diviTest, Special) { testSpecial(&func); }                \
+  TEST_F(LlvmLibc##Name##diviTest, Special) { testSpecial(&func); }            \
   static_assert(true, "Require semicolon.")

Copy link
Contributor

@PiJoules PiJoules left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we're almost there. Just a couple more changes to test I'd like to see for more confidence.

EXPECT_EQ(func(-32768,32768), FRACT_MIN);
EXPECT_EQ(func(32767,32768), FRACT_MAX);
EXPECT_EQ(func(INT_MAX,INT_MAX), 1.0r);
EXPECT_EQ(func(INT_MAX-1,INT_MAX), 0.99999999r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since 0.99999999r also can't fit precisely, let's use the EXPECT_LT epsilon approach used above.

Comment on lines +49 to +50
EXPECT_EQ(func(-32768,32768), FRACT_MIN);
EXPECT_EQ(func(32767,32768), FRACT_MAX);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this will only be true if the number of fractional bits for a fract type is 16, but it can vary depending on the platform, so let's maybe instead derive these number of fractional bits. Something like

EXPECT_EQ(func(-1 << FRACT_FBIT, 1 << FRACT_FBIT), FRACT_MIN);
EXPECT_EQ(func((1 << FRACT_FBIT) - 1, 1 << FRACT_FBIT), FRACT_MAX);

/* Expecting 0 here as fract is not precise enough to
* handle 1/INT_MAX
*/
EXPECT_EQ(func(1, INT_MAX), 0.r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to before, rounding can be either up or down, so I think we also want to use the epsilon approach here.

EXPECT_EQ(func(0,-10), 0.r);
EXPECT_EQ(func(-32768,32768), FRACT_MIN);
EXPECT_EQ(func(32767,32768), FRACT_MAX);
EXPECT_EQ(func(INT_MAX,INT_MAX), 1.0r);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is correct but it had me second guessing a little because of some nuances with the spec I had to look up.

Per 7.18a.6.1

For functions returning a fixed-point value, the return value is saturated on overflow. 

so rdivi(INT_MAX,INT_MAX) should definitely return FRACT_MAX, and per Clause 6.4.4:

The value of a constant shall be in the range of representable values for its type, with exception for constants of a fract type with a value of exactly 1; such a constant shall denote the maximal value for the type.

so 1.0r should also resolve to FRACT_MAX hence the comparison is correct. Just to make things clearer, could you add these as comments above this check, or alternatively make the RHS FRACT_MAX while keeping the rdivi saturation comment?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tangentially, it also might be good to add some tests which would trigger the saturated overflow case. So maybe some integral values which would resolve to either less than zero or greater than or equal to one which would saturate to zero and FRACT_MAX.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants